﻿
using System;
using System.Collections.Generic;
using System.Linq;

namespace LxCFD
{
    public class SIMPLEBase : LinearSolver, ISolvableUnit// SolvableUnitBase
    {
        public List<VelocityVector> Vels = new List<VelocityVector>();
        public List<NSEquation> Us = new List<NSEquation>();
        public List<NSEquation> Vs = new List<NSEquation>();
        public List<NSEquation> Ws = new List<NSEquation>();
        public List<PressureCorrectionBasic> PCors = new List<PressureCorrectionBasic>();
        public List<PressureSIMPLE> Ps = new List<PressureSIMPLE>();
        public List<ConnectBase> Connects = new List<ConnectBase>();//暂时不用，占位

        protected volatile bool isExistPressureReference;
        protected double PressureReferenceCorrection;
        protected int CenterZoneIndex;

        public void CreateUnit(List<PressureCorrectionBasic> scalars, List<ConnectBase> connects,SolverTypes solverType)
        {
            this.PCors = scalars;
            foreach (var pc in PCors)
            {
                Vels.Add(pc.VelocityVector);
                Us.Add(pc.VelocityVector.U);
                Vs.Add(pc.VelocityVector.V);
                Ws.Add(pc.VelocityVector.W);
                Ps.Add(pc.PressureSIMPLE);
            }
            this.Connects = connects;
            SolverType = solverType;
            SetAMGSolver(scalars, connects);
            setReferencePressure();
        }

        public void setReferencePressure()
        {
            isExistPressureReference = false;
            foreach (var v in Vels)
            {
                if (v.PressureInOuts.Count > 0)
                {
                    isExistPressureReference = true;
                }
            }
            if (!isExistPressureReference)
            {
                Vector3d centerCal = new Vector3d();
                double totalVolumn = 0;
                foreach (var p in Ps)
                {
                    centerCal += p.Zone.Volumn * p.Center;
                    totalVolumn += p.Zone.Volumn;
                }
                centerCal.MultiplyToSelf(1 / totalVolumn);
                double minlenthSquare = double.MaxValue;
                for (int i1 = 0; i1 < Ps.Count; i1++)
                {
                    double ls = (Ps[i1].Center - centerCal).CalLengthSquared();
                    if (ls < minlenthSquare)
                    {
                        minlenthSquare = ls;
                        CenterZoneIndex = i1;
                    }
                }
            }
        }

        protected void updateMomentMatrix()
        {
            //1.求解动量方程
            MultiThreads.RunCalTasks((taskIndex) =>
            {
                foreach (var v in Vels) { v.SetUVW_x(taskIndex); }
            });
            MultiThreads.RunCalTasks((taskIndex) =>
            {
                foreach (var v in Vels)
                {
                    v.U.CalGradient(taskIndex);
                    v.V.CalGradient(taskIndex);
                    v.W.CalGradient(taskIndex);
                }
            });
            MultiThreads.RunCalTasks((taskIndex) =>
            {
                foreach (var v in Vels)
                {
                    v.U.SetAb_Step1(taskIndex); v.V.SetAb_Step1(taskIndex); v.W.SetAb_Step1(taskIndex);
                    v.AddViscositySource_LoopFace(taskIndex);
                }
            });
            MultiThreads.RunCalTasks((taskIndex) =>
            {
                foreach (var v in Vels)
                {
                    v.U.SetAb_Step2(taskIndex); v.V.SetAb_Step2(taskIndex); v.W.SetAb_Step2(taskIndex);
                    v.AddSource_LoopCell(taskIndex);//每个线程的区域是固定的，不会相互干涉
                }
            });
        }
        protected void updatePCors()
        {
            //2.Rhie-Chow求解massImbalance
            MultiThreads.RunCalTasks((taskIndex) =>
            {
                foreach (var v in Vels) { v.CalFaceMassFlux(taskIndex); }
            });
            MultiThreads.RunCalTasks((taskIndex) =>
            {
                foreach (var v in Vels) { v.CalMassImbalance(taskIndex); }
            });
            //3.求解压力修正
            MultiThreads.RunCalTasks((taskIndex) =>
            {
                foreach (var pc in PCors) { pc.Set_x_Bound(taskIndex); }
            });
            MultiThreads.RunCalTasks((taskIndex) =>
            {
                foreach (var pc in PCors) { pc.SetAb_Step1(taskIndex); }
            });
            MultiThreads.RunCalTasks((taskIndex) =>
            {
                foreach (var pc in PCors) { pc.SetAb_Step2(taskIndex); }
            });
        }
        protected void updatePressure()
        {
            //5.压力场,流场更新
            MultiThreads.RunCalTasks((taskIndex) =>
            {
                foreach (var p in Ps) { p.Predict_bound(taskIndex); }
            });
            if (!isExistPressureReference)
            {
                MultiThreads.RunCalTasks((taskIndex) =>
                {
                    foreach (var p in Ps) { p.CalPressureArea(taskIndex); }
                });
                {
                    double iniPA = 0, boundArea = 0;
                    foreach (var z in Ps)
                    {
                        iniPA += z.InitiialPressureProductBoundArea;
                        boundArea += z.BoundArea;
                    }
                    ThreadsMath.GetSum(out double currentPressureProductArea, out double sum2, out double sum3);
                    PressureReferenceCorrection = (iniPA - currentPressureProductArea) / boundArea;
                }
                MultiThreads.RunCalTasks((taskIndex) =>
                {
                    foreach (var p in Ps) { p.CorrectPressureByReferenct(PressureReferenceCorrection, taskIndex); }
                });
            }
            MultiThreads.RunCalTasks((taskIndex) =>
            {
                foreach (var p in Ps) { p.CalGradient(taskIndex); }
            });
            MultiThreads.RunCalTasks((taskIndex) =>
            {
                foreach (var v in Vels) {
                    v.CalResidual(taskIndex);
                    v.U.CalResidual(taskIndex);
                    v.V.CalResidual(taskIndex);
                    v.W.CalResidual(taskIndex);
                }
            });
            MultiThreads.RunCalTasks((taskIndex) =>
            {
                foreach (var v in Vels) { v.UpdateUVWBound(taskIndex); }
            });
        }
        public virtual void Solve_Steady()
        {
            updateMomentMatrix();
            MultiThreads.RunCalTasks((taskIndex) =>
            {
                foreach (var pc in PCors) { pc.CalDc_ADiagNotRelaxed(taskIndex); }
            });
            RelaxAndSolve(Us, Connects);
            RelaxAndSolve(Vs, Connects);
            RelaxAndSolve(Ws, Connects);
            updatePCors();
            solvePCors();
            //4.进行修正
            MultiThreads.RunCalTasks((taskIndex) =>
            {
                foreach (var pc in PCors) { pc.ApplyPCorToUVWP(taskIndex); }
            });
            MultiThreads.RunCalTasks((taskIndex) =>
            {
                foreach (var v in Vels) { v.CalMassImbalance(taskIndex); }
            });
            updatePressure();
            //其他场
        }
        public virtual void Solve_ImplicitTransient()
        {
            updateMomentMatrix();
            MultiThreads.RunCalTasks((taskIndex) =>
            {
                foreach (var v in Vels)
                {
                    v.U.TransData.SetAb_ImplicitTrans(taskIndex);
                    v.V.TransData.SetAb_ImplicitTrans(taskIndex);
                    v.W.TransData.SetAb_ImplicitTrans(taskIndex);
                }
            });
            MultiThreads.RunCalTasks((taskIndex) =>
            {
                foreach (var pc in PCors) { pc.CalDc_ADiagNotRelaxed(taskIndex); }
            });
            RelaxAndSolve(Us, Connects);
            RelaxAndSolve(Vs, Connects);
            RelaxAndSolve(Ws, Connects);
            updatePCors();
            solvePCors();
            MultiThreads.RunCalTasks((taskIndex) =>
            {
                foreach (var pc in PCors) { pc.ApplyPCorToUVWP(taskIndex); }
            });
            MultiThreads.RunCalTasks((taskIndex) =>
            {
                foreach (var v in Vels) { v.CalMassImbalance(taskIndex); }
            });
            updatePressure();
        }
        protected virtual void solvePCors()
        {
            if (!isExistPressureReference)
            {
                PCors[CenterZoneIndex].SetAb_PCor0Point(Ps[CenterZoneIndex].CenterCellIndex);
            }
            RelaxAndSolve(PCors, Connects);
            MultiThreads.RunCalTasks((taskIndex) =>
            {
                foreach (var pc in PCors) { pc.Update_x_Bound(taskIndex); }
            });
            MultiThreads.RunCalTasks((taskIndex) =>
            {
                foreach (var pc in PCors) { pc.CalGradient(taskIndex); }
            });
        }
        public void SaveOldValue()
        {
            foreach (var v in Vels)
            {
                v.U.TransData.PrepareToSave();
                v.V.TransData.PrepareToSave();
                v.W.TransData.PrepareToSave();
                v.TransVelData.PrepareToSave();
            }
            MultiThreads.RunCalTasks((taskIndex) =>
            {
                foreach (var v in Vels)
                {
                    v.U.TransData.SaveOldValue(taskIndex);
                    v.V.TransData.SaveOldValue(taskIndex);
                    v.W.TransData.SaveOldValue(taskIndex);
                    v.TransVelData.SaveOldValue(taskIndex);
                }
            });
        }

        public void Solve_ExplicitTransient()
        {
            throw new NotImplementedException();
        }
    }
}


//public override List<SIMPLEBase> SplitByConnectivity(ConnectivityAnalysisTool connectivityAnalysis)
//{
//    List<SIMPLEBase> Units = new List<SIMPLEBase>();
//    List<ConnectivityAnalysisTool.Range> pools = connectivityAnalysis.SplitByConnectivity(Vels, connects);
//    foreach (var p in pools)
//    {
//        SIMPLEBase u;//= new SIMPLEBaseUnit();
//        if (PCors[p.zonesIndex[0]].solvingMethods == SolvingMethods.PISO)
//        {
//            u = new PISOUnit();
//        }
//        else
//        {
//            u = new SIMPLE();
//        }
//        foreach (int z in p.zonesIndex)
//        {
//            u.AddVelocityAndPressure(Vels[z], PCors[z], Ps[z]);
//        }
//        foreach (int c in p.contactsIndex)
//        {
//            //u.AddConnect(connectTs[c]);
//        }
//        //u.addAMGSolver(TempStor.Vels, TempStor.connects);
//        //u.SolverType = this.solverType;
//        //u.setReferencePressure();
//        Units.Add(u);
//    }
//    return Units;
//}
//public override void InitializeLineSolver(SolverTypes solverType)
//{
//    base.InitializeLineSolver(solverType);
//    addAMGSolver(Vels, connects);
//    setReferencePressure();
//}
// public SIMPLEBaseUnit TempStor { get; set; }
//public override void Clear()
//{
//    base.Clear();
//    TempStor = new SIMPLEBaseUnit();
//}
//public void AddVelocityAndPressure(VelocityVector vel, PCorBase PCor, PressureInSimple P)
//{
//    Vels.Add(vel);
//    Us.Add(vel.U);
//    Vs.Add(vel.V);
//    Ws.Add(vel.W);
//    PCors.Add(PCor);
//    Ps.Add(P);
//}